/*******************************************************************************

                                APODFC

  Title:        Apodize a file of fcalc data and prepare Wilson-style plots.

  Author:       Hanna Szoke

  Date:         10/05/00: Minor changes to bring into line with apodfo.c.
  		04/10/03: Good for both gridded and list Fc files.

*******************************************************************************/

#include "util.h"
#include "cellparams.h"
#include "dims.h"	
#include "apodize.h"
#include "sp.h"

static	int	Ntot ;			/* # masked-in structure factors */
static  real	*Int ;			/* full set of (amp^2)	*/
static  real	*ksq ;			/* full set of (1/d^2)  */
static  real	*sigsqI ;

static	int	apodfc_startup() ;
static	void	apply_apod_fcList() ;	/* use the delta that was calculated  */
static	void	prepare_fcList() ;		/* make weighted R^2 */

void
apodfc_main(int argc, char *argv[])
{
	/* in apodutil.c:       */

     void	allocate_Bin_arrays(int) ;
     void	do_binning(int, real *, real *, real *) ;	
     real	get_delta(real) ;
     void	init_apod_arrays(int, real *, real *, real *) ;
     real	linfit(int, real *, real *) ;		
     void	make_wilson(real *) ;
     void	make_selected_wilson(char *) ;
     void	show_wilson_plots() ;

     void	write_HDFfcalc(char *) ;
     void       writeFclist(char *) ;

     char	out_filename[MAXSTRING] ; /* name of output file for x/y data */
     real	delta1 ;
     real	dstdev1 ;
     real	B1 ;

  /*******************************************************
  Check for right ballpark; identify Fcalc file name.
  *******************************************************/

     if (argc > optind+3) 	/* too many arguments */
	  ballpark(caller) ; 
     else if (argc == optind+3)
          sprintf(sf_filename, "%s", argv[optind+2]) ;
     else {
	  prompt("\nPlease enter the name of an fcalc file: ") ;
          sscanf(terminp, "%s", sf_filename) ;
	  sprintf(message, "Running %s on fcalc file %s", caller, sf_filename) ;
	  printTwice(message) ;
     }

     strip_path(sf_filename, pwd_name) ;

     hello(caller) ;

  /*******************************************************
  Read the input fcalc file using standard read function.
  Set up arrays in (h,k,l) reciprocal lattice triplets,
  Prepare a list of non-zero intensities, dummy sigma
  squared*intensities, and nusq's.  
  *******************************************************/

     Nbin = apodfc_startup() ;

     init_apod_arrays(Ntot, Int, ksq, sigsqI) ;
     prepare_fcList() ;

  /*******************************************************
  Bin the information, gathering logs of average binned 
  intensities & weights, 
  *******************************************************/

     allocate_Bin_arrays(Nbin) ;
     do_binning(Ntot, Int, ksq, sigsqI) ;

  /*******************************************************
  Do a linear fit to establish slope and intercept.
  Prepare for display the original binned log(I^2) and
  the linearly fitted version.
  *******************************************************/

     B1 = linfit(Nbin, w_lnIsq, apodwt) ; 
     delta1 = get_delta(B1) ; 
     make_wilson(&dstdev1) ;

     show_wilson_plots() ;
     delta = delta1 ;

  /*******************************************************
  Replace fc input file by apodized version, insofar as 
  delta > 0.
  *******************************************************/

     if (delta > 0) {
          apply_apod_fcList() ;
          if (HDFOUT){
               extend_filename(pwd_name, out_filename, "_apo");
               write_HDFfcalc(out_filename) ;
          }
          else{
               extend_filename(pwd_name, out_filename, "_apo.fcalc");
               writeFclist(out_filename) ;
          }
          init_apod_arrays(Ntot, Int, ksq, sigsqI) ;
          prepare_fcList() ;

          do_binning(Ntot, Int, ksq, sigsqI) ;
     }
     sprintf(out_filename, "%s_wil", pwd_name) ;
     make_selected_wilson(out_filename) ;

     printTwice("\nPLEASE NOTE:") ;
     sprintf(message, " The Wilson file '%s' should be used for scaling your fobs.",
     out_filename) ;
     printTwice(message) ;
}

void prepare_fcList()
{
     real	fc_nusq(int) ;
     real	ampsq ;
     real	Lreal, Limag ;
     int	n ;

     for (n = 1; n < Nfc; n++) {
	     
          Lreal = *(Fc_relist + n) ;
	  Limag = *(Fc_imlist + n) ;

	  ampsq = Lreal*Lreal + Limag*Limag ;

          *(Int+n) = ampsq ;
	  *(ksq+n) = fc_nusq(n) ;
          *(sigsqI+n) = 4. * ampsq ; 
     }			

     return ;
}


void	apply_apod_fcList() 	
{	  
     real	fc_nusq(int) ;
     COMPLEX	oneitem ;
     real	amp ;
     real	phase ;
     int	n ;

     for (n = 1; n < Nfc; n++) {
	     
          oneitem.re = *(Fc_relist + n) ;
	  oneitem.im = *(Fc_imlist + n) ;

          getAmpPhase(&oneitem, amp0, &amp, &phase) ;

          amp *= exp(-(delta * delfac * fc_nusq(n))) ;

          *(Fc_relist + n) = amp * cos(DTOR*phase) ;
          *(Fc_imlist + n) = amp * sin(DTOR*phase) ;
     }			

     return ;
}

  /*******************************************************
  Read the input fcalc file using standard read function.
  Set up arrays in (h,k,l) reciprocal lattice triplets,
  *******************************************************/

int apodfc_startup()
{
     void	a_readInput() ;		/* read input parameters */
     void	alloc_fc() ;
     int	HDFdatalength(char *, char *) ;
     int	read_HDFfcalc(char *) ;
     int	readFcList(char *) ;
     void	fsetDsqLimit(int, char) ;
     int	n ;

     a_readInput() ;

     strcpy(message, "apodfc_startup") ;

     if (HDFIN) {
          Nfc = HDFdatalength(sf_filename, "/fcalc/Nfc") ;
          alloc_fc(Nfc) ;
          Ntot = read_HDFfcalc(sf_filename) ;
     }
     else {
          Nfc = filelength(sf_filename) ;
          alloc_fc(Nfc) ;
          Ntot = readFcList(sf_filename) ;
     }

     if (Ntot == 0) 
	  Fatal("Cannot apodize an empty file!") ;

     Int = (real *) e_malloc(Ntot*sizeof(real), message) ; 
     sigsqI = (real *) e_malloc(Ntot*sizeof(real), message) ; 
     ksq   = (real *) e_malloc(Ntot*sizeof(real), message) ; 

     fsetDsqLimit(Nfc, FALSE) ;

     n = (dsq_limit + 0.5*binwidth)/ binwidth + 1 ;

     sprintf(message, "# of bins set to %d, average bin size is %g.", 
		n, binwidth) ;
     printTwice(message) ;

     return (n) ;
}
